function [LS, LD, rho_grid]=LSLD(rho_grid,func_var,pi_m,Nb,y)
rq_ratio = func_var.rq_ratio;cost=func_var.cost;
beta = func_var.beta;df_inv = func_var.df_inv;
xi_func=@(x)(1-rq_ratio).*max(1+x,1/(1+pi_m))+rq_ratio/(1+pi_m)-cost;
D_bar = max(y(:,1));
psi_func=@(x)interp1(y(:,1),y(:,3),min(x,D_bar)).*(x<=D_bar)+beta.*(x>D_bar);
dpsi_func=@(x)interp1(y(:,1),y(:,4),min(x,D_bar)).*(x<=D_bar);

eq_D_func = @(x,y)xi_func(y).*(dpsi_func(x).*x/Nb+psi_func(x))-1;
if rho_grid(1)<1/(1+pi_m)-1
    ID = find(rho_grid>1/(1+pi_m)-1,1,'first');
    rho_grid = [rho_grid(1:ID-1),1/(1+pi_m)-1,rho_grid(ID:end)];
end

for i=1:length(rho_grid)
    [D_grid(i),fval(i)] = fzero(@(x)eq_D_func(x,rho_grid(i)),[1e-4,D_bar]);
  
    LS(i) = (1-rq_ratio)*D_grid(i)*psi_func(D_grid(i));
end

if rho_grid(1)<1/(1+pi_m)-1
    ID = find(rho_grid<1/(1+pi_m)-1,1,'last');
    LS(1:ID)=0;
end

LD = df_inv(rho_grid);

% %debug
% z_grid = linspace(0.5,1,100);
% i=91
% figure
% hold on
% plot(z_grid,eq_D_func(z_grid,rho_grid(i)))
% plot(z_grid,eq_D_func(z_grid,rho_grid(i+1)))
% plot(z_grid, z_grid*0)
% 

